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This work identifies geometric effects on dynamics due to nonadiabatic couplings in Born Op- 
penheimer systems and provides a systematic method for deriving corrections to mixed quantum- 
classical methods. Specifically, an exact path integral formulation of the quantum nonadiabatic 
dynamics of Born Oppenheimer systems is described. Stationary phase approximations to the prop- 
agator for full quantum dynamics are derived. It is shown that quantum corrections to mixed 
quantum classical methods can be obtained through stationary phase approximations to the full 
quantum dynamics. A rigorous description of the quantum corrections due to electronic nonadi- 
abatic coupling on the nuclear dynamics within the Ehrenfest framework is obtained. The fewest 
switches surface hopping method is shown to be obtained as a quasiclassical approximation to the 
dynamics and natural semiclassical extensions to include classically forbidden nonadiabatic transi- 
tions are suggested. 

I. Introduction 

The Born-Oppenheimer approximation^ is fundamental to studies of a wide variety of molecular systems. How- 
ever, violations of the Born-Oppenheimer approximation^ are ubiquitous in many molecular and condensed phase 
phenomena, and are responsible for several relaxation processes. Despite their importance, fully quantum treatments 
of nonadiabatic corrections to Born Oppenheimer dynamics are limited and no perturbative formulation exists to 
account for them. 

It is practically impossible to compute the full quantum nonadiabatic dynamics for realistic systems due to the 
prohibitive computational cost associated with such a calculation. To circumvent this difficulty, a class of mixed 
quantum classical methods have been invented^iii^i^ili^i^iiS. These methods can be broadly classified into mean- 
field^i2i»i^ and trajectory surface hopping methods^i^. Mean field methods originate from work by Ehrenfest'^, where 
the nuclear degrees of freedom, treated as classical mechanical variables, experience a potential that is the expectation 
value of the electronic Hamiltonian with respect to the instantaneous electronic wavefunctions obtained by solving 
the electronic Schrodinger equation. On the other hand, surface hopping methods^i^ treat the nuclear dynamics 
as being localized to a single Born-Oppenheimer surface at any time, with nonadiabatic couplings between surfaces 
leading to instantaneous jumps, or "hops" between Born-Oppenheimer surfaces. Common to all mixed quantum- 
classical methods is the identification of the nuclear (or slow) degrees of freedom as behaving as essentially classical 
quantities, while interacting with a quantum system consisting of the electronic degrees of freedom. This identification 
is usually justified by a mixture of mathematical and phenomenological arguments. As a result, while these methods 
have been investigated by intuitive tests and computational simulations over the years, there is no rigorous theory 
to understand situations where this separation into classical and quantum subsystems breaks down and quantum 
effects on the classical degrees of freedom have to be considered. This problem is also evident in studies of the 
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dynamics of molecular systems on higher-dimensional Born-Oppenheimer surfaces, where geometric features, like 
conical intersections between surfaces, have a significant physical effect on the dynamics. Thus, given the somewhat 
ad-hoc nature of most mixed quantum classical methods, a central problem is to develop a framework which would 
incorporate these methods as an approximation and would enable the systematic study of corrections to the mixed 
quantum classical picture of nonadiabatic dynamics. 

This work aims to systematically construct mixed quantum classical methods as scmiclassical limits of quantum 
nonadiabatic dynamics. Such an approach enables the systematic improvement of mixed quantum classical meth- 
ods to include further quantum corrections as well as the development of new and hybrid mixed quantum classical 
methods. Exact path integral descriptions of the quantum nonadiabatic dynamics are provided and stationary phase 
approximations for the dynamics are derived. The resulting dynamical picture is a consistent means of obtaining a 
mixed quantum classical approximation to the nonadiabatic dynamics. It is demonstrated that both the Ehrenfest ap- 
proach and the surface hopping approach correspond to stationary phase approximations of the quantum nonadiabatic 
dynamics, when studied in complementary representations of the electronic system. Furthermore, the approximate 
equations of motion obtained in the Ehrenfest picture contain additional contributions to the effective nuclear force. 
These contributions are a direct, and nontrivial consequence of the geometry of the Born-Oppenheimer potential 
energy surfaces. For systems with only one nuclear dimension, these contributions are zero, and the dynamics is 
of pure Ehrenfest type, while for higher dimensional surfaces, they are non-zero, and represent corrections to the 
quasiclassical dynamics due to the topological features of the surface, like the presence of conical intersections and 
degeneracies. In the complementary picture from which surface hopping methods are obtained as an approximation, 
quantum corrections to the nonadiabatic nuclear dynamics can be built in by considering scmiclassical expansions of 
the nuclear propagator in the Born-Oppenheimer approximation. In addition to this, the approach outlined in this 
paper provides a new interpretation of nonadiabatic effects as a consequence of the geometry of the electronic Hilbert 
space, and thus as arising due to a generalization of Berry's phase for parametrized quantmn systems. 

The outline of the paper is as follows. In Sec.II, a path integral description of the quantum nonadiabatic dynamics 
is introduced. Issues concerning path integral approaches to nonadiabatic dynamics are discussed and the derivation 
of the approach used in this paper are outlined in Sec. II. A. In Sec. II. B, the path integral description is derived. It is 
based on integrating out the electronic degrees of freedom using a coherent state basis. Furthermore, the path integral 
so obtained is compared to other path integral formulations constructed within a Born-Oppenheimer basis. It is shown 
in this section that nonadiabatic couplings between states are due to the geometry of the electronic Hilbert space, and 
can be regarded as generalizations of Berry's phase. In Sec. Ill, the path integral approach is analysed and stationary 
phase approximations are derived. Sec. III. A demonstrates that the stationary phase approximation to the coherent 
state path integral corresponds to the Ehrenfest formulation of mixed quantum classical dynamics. Furthermore, 
geometric corrections to the nuclear force in the Ehrenfest approach are obtained. These corrections are related 
to the phase coherence between different states, and are most significant when different Born Oppenheimer energy 
surfaces become close in energy and are strongly coupled due to nonadiabatic effects. Sec. III. B studies stationary 
phase approximations to the complementary path integral description constructed using a Born-Oppenheimer basis 
for the electronic system. It is shown that the equations of motion so obtained correspond in the scmiclassical limit, 
to the fewest switches surface hopping approach to nonadiabatic dynamics. The theory in Sec.III.B also points to the 
construction of scmiclassical propagators that can extend descriptions of nonadiabatic dynamics beyond the mixed 
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quantum-classical limit. Sec. III. C describes semiclassical quantization rules to approximate the eigenvalue spectrum 
for the full quantum dynamics. The paper then concludes with a discussion of the results and future work in Sec. IV. 

II. Theory 

The problems with path integral approaches to quantum nonadiabatic dynamics are briefly discussed, and a short 
outline of the consequent derivation of the path integral based theory is given. The detailed mathematical treatment 
of a path integral approach based on coherent states is then subsequently discussed. 

It should also be noted that this text uses the units convention where Planck's constant is set to one, i.e h — 1 in 
this paper. 

A. Preliminaries 

The nonadiabatic corrections to the Born Oppenheimer approximation can be viewed as a result of a coupling 
between a set of "fast" or bath degrees of freedom and slow or "system" degrees of freedom. However, the mechanisms 
that result from this coupling are different from the usual problem of a system interacting with a dissipative bath, 
since the Hilbert space enclosing the bath degrees of freedom in the Born Oppenheimer theory is parametrized by the 
slow degrees of freedom. As a consequence, the time dependent evolution of such a system involves complex geometric 
phase effect o^-'^i-'^^'^'^i^^'^^i^^ . Pechukas^i^ first developed a path integral prescription for nonadiabatic dynamics and 
derived a set of nonlocal equations of motion for the semiclassical nuclear path. Later wor k^^i^^ on a path integral 
description clarified the interpretation of the nonadiabatic coupling vector potentials as nonabelian gauge fields. 
However, these descriptions are not easily amenable to approximations which can be studied computationally. In the 
Pechukas formulation, the effect of the electronic degrees of freedom on system dynamics enters through a nonlocal 
potential which is a function of the paths taken by the system, and similarly in the formulation of Moody, Shapere and 
Wilczeki^iii, nonadiabatic couplings between Born-Oppenheimer states enter as path ordered matrix exponentials. 
Thus, an approach that adopts the path integral framework for the exact nonadiabatic dynamics, but which is 
also amenable to systematic and computationally tractable approximations would be of considerable theoretical and 
practical utility. 

The problems with regards to nonlocality and path ordered exponentials encountered in a path integral description 
of the quantum nonadiabatic dynamics can be partially resolved by an appropriate choice of the path integral. To 
this end, a local basis of continuous, A^-particle coherent state Slater determinants is chosen to describe the electronic 
system. The use of coherent states enables the straightforward construction of stationary phase approximations to 
the quantum dynamics and corresponding semiclassical quantization rules. Furthermore, since the propagator can be 
written as a functional integral over a continuous coherent state space, issues relating to path ordering and nonlocality 
become considerably simplified. 

The coherent state basis is parametrized by the nuclear coordinates and is overcomplete. Since the basis is 
parametrized by the nuclear coordinates, a local Hilbert space is defined at each point on a given nuclear path. 
Thus, intuitively, motion along a nuclear path involves the following: The first contribution is due to the propagation 
of the nuclear part of the Hamiltonian, and the second is due to a rotation of the Hilbert space from one point into 
the next one on the path (Fig 1). This rotation takes basis states at one point into the basis states at the next, and 
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hence changes the electronic contribution to the effective nuclear Hamiltonian. In addition, the rotation of the basis 
in general involves "twists" and skews resulting in torsional contributions to the force acting on the nuclear degrees 
of freedom. 

To explicitly describe this, the derivation of the coherent state path integral proceeds as follows. The matrix 
elements of the time propagation operator for the full (nuclear plus electronic) Hamiltonian between initial and final 
nuclear positions are written down as a path integral. The time for propagation is broken up into a large number of 
(infinitesimal) time segments and the propagation operator is written in a Trotter product expansion. Complete sets 
of nuclear position eigenstates are then introduced for each time segment, and an infinitesimal propagator between 
successive nuclear positions is written down. The resulting expression is a path integral of an operator acting on the 
electronic degrees of freedom over each nuclear path. For a given nuclear path, this quantity simply corresponds to 
the electronic propagation operator parametrized by the nuclear path. 

The matrix elements of the quantum nonadiabatic propagation operator between initial and final electronic states 
and nuclear positions are then evaluated through the matrix elements of this path dependent quantity between initial 
and final electronic states. This is done by adding overcomplete sets of electronic A'' particle coherent state Slater 
determinants at each time slice and evaluating matrix elements of the infinitesimal propagation operator. The overlap 
between these basis states belonging to consecutive timesteps gives rise to the nonadiabatic coupling contribution to 
the dynamics. The resulting path integral is then further analysed. 



B. Derivation 



The full Hamiltonian for a system with nuclear and electronic degrees of freedom can be written as 

F = if„(R)+ife(r,R) (1) 

The nuclear Hamiltonian ff„(R) is given by 

^"W = E^ + f^W (2) 

1=1 ^ 

The coordinates R = (.Ri, ...Rk) denote the K nuclear degrees of freedom with momenta given by P = (Pi,P2, ...P^)- 
The potential energy function, U (R) is the bare nuclear potential. The electronic Hamiltonian can be similarly defined 
and the electronic coordinates, {fi,r2...rN}, are elements of the vector r. The electronic potential energy y(r, R) is 
parametrically dependent on the nuclear positions. 



N 
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ffe(r,R) = ^^ + F(r,R) (3) 

For simplicity of presentation, the nuclei and electrons are assumed to be spinless. This assumption is reasonable for 
the nonadiabatic effects considered here, and the treatment that follows can be extended to include spin effects. Also 
it is assumed that the nuclear masses, M/ are all equal and have the value M. This is solely for notational convenience 
and does not significantly affect the conclusions of this work. 

A Born Oppenheimer basis of eigenstates and eigenvalues is defined by 

[ffe(r,R)] I M/^;R) = E^{K) I M/^;R) (4) 
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The unitary time propagation operator for this Hamiltonian can be expanded in a Trotter product expansion as 
described below. 

P - - - T>2 



n '^^p 



k=l 



exp 



^2M 



exp 



U) 



(5) 



The propagation operator for a given set of initial and final nuclear positions Q and Q' can be written as 
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(6) 



where Ri = Q and Rp = Q'. It is important to note that the matrix elements of the propagation operator described 
in Eq.(6) are still operators acting on the electronic Hilbert space. This is because the matrix elements are a partial 
trace with respect to the eigenstates | Q') and | Q) of the nuclear position operator R. This partial trace can be 
rewritten as a path integral over the set of nuclear paths R(t). 



(Q' I e-'^* I = / ^R(^)e^'^"f^'"'lG„a[R(r)] 



(7) 



S'„[R(i)] is the bare nuclear action corresponding to the potential [/(R). The quantity Gna is an operator quantity in 
the electronic degrees of freedom, and is a nonlocal functional of the nuclear paths. Physically, the operator G„a[R-(''')] 
includes the "connection" between electronic Hilbert spaces defined at each point on a given nuclear path R(t). Thus, 
a complete solution of the nonadiabatic problem would require that electronic matrix elements of Gna be evaluated. 
To perform such an evaluation, overcomplete sets of antisymmetric N electron Slater determinant coherent state 
wavefunctions are introduced at each nuclear position along a given nuclear pat h^^i^^ . 

A general N electron antisymmetric wavefunction can be written as a linear combination of Slater determinants. 
This basis of N electron Slater determinants, | <j}) =\ 01, ■■■4>n) is overcomplete and has the closure property: 

H0P 



1 = 



(8) 



Here, a local basis of TV electron coherent state Slater determinant a-*^^'^^ is defined at each point on a given nuclear 
path. Hence, this local electronic basis is parametrized by the nuclear position to account for the fact that a local 
electronic Hilbert space exists at each nuclear position. For such a basis defined for a nuclear position R, the 
completeness relation reads: 

*(R)d0(R) 



1 



(R) 



(R))(0(R) 



For notational convenience, the integration measure in the overcompleteness relation will be written as 



*(R)d(/)(R) 



(R.)r 



with the differentials and the absolute value of (?I>(R) defined by 

N 



(9) 



(10) 



N 
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The parametrization of the electronic Hilbert space by the nuclear position is responsible for geometric phase effects. 
As will be shown below, these effects enter through the definition of a local basis of wavefunctions at each nuclear 
position. 
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The operator Gna can be explicitly defined to be 

p 

Gnamt)] = n exp [-ieH,{r, RJ] (12) 
fe=i 

This operator is a path ordered quantity because the electron Hamiltonian does not commute with itself when 
evaluated at different points on a nuclear path. To evaluate the matrix elements of this operator, overcomplete sets 
of A/'-electron Slater determinant wavefunctions are introduced at each timestep. On doing so, the matrix elements 
of Gna between an initial electronic state | Q) and a final state | Q') can be written as 

r ^ 

($°;Q' I Gna I $^;Q> = / l[dn[^l,jiRj)]dn[ct>jiRj)]{<P°;Q' \ V(Q')>(V'(Q') I e-^=(-'Q > | <^(Q')> 

(.^(QO I Vi(Ri))(V'i(Ri) I e--^^(--'«^) I <^i(Ri))...(.^p(Rp) I $/;Q) 

(13) 

Two sets of overcomplete basis states are introduced at each timestep, because the time evolution operator is in 
general not diagonal in this basis of states, and the overlap between basis states corresponding to different positions 
on the nuclear path is not unity. More specifically, the coherent states inserted at successive timestcps, tk,tk+i, are 
parametrized by different sets of nuclear positions, Rfe,Rfe-|-i. Since there is no constraint on the basis states at 
different nuclear positions, all possible matrix elements of the time evolution operator between these basis states will 
have to be considered. Furthermore, to evaluate this off-diagonal matrix element, an additional set of basis states 
parametrized by the position R^ is introduced at time tk+i- Since the path integral is discontinuous in the coherent 
state space, all possible overlaps of this new basis element with the basis element parametrized by Rfe+i will have 
to be considered, and the matrix elements of the time evolution operator are not constrained to be diagonal. This 
lack of constraint is unlike the cases encountered with more traditional coherent state path integrals where only 
diagonal matrix elements between coherent states suffice, and is a consequence of the parametrization by the nuclear 
coordinates. 

It is clear that the above path integral involves matrix elements of the form 

A[V'(Rfc), <^(Rfe)] = (V(Rfc) I exp {-ieH{v, RJ) | 0(Rfc)) (14) 

If the exponent on the right hand side of Eq.(14) is expanded in powers of e and terms higher than the first order are 
neglected, the matrix elements A can be evaluated to be approximately 

A[V'(Rfe), <^(Rfe)] = (V'(Rfc) I [l - ieH{v, Rfc)] | cj,{Kk)) (15) 

Making the same assumption as before, the right hand side of Eq.(15) can be reexponentiated to yield 

A[V-(R.), <^(R.)] = {^k) I <^(R.)) exp [ - (16) 

In addition to the matrix elements A, the path integral contains "connection" factors of the form 

Tfc = (<^(Rfc_i) I V(Rfc)) (17) 



These connection factors are geometric in character and arc due to the parallel transport of the local electronic basis 
from one nuclear position to the next. Thus, these factors correspond to a generalization of the canonical, or Berry's 
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phase for simple, parameter dependent Hamiltonians^'^. To obtain a more familiar form for this phase dependence, 
the connection factors can be rewritten as follows: 



The difference (5(/)(Rfe) is the difference between the wavefunctions cj) evaluated at successive time steps. 



S^Rk 



(Rfe)-0(Rfe,i) 



(18) 



(19) 



On exponentiating and making the assumption that the basis wavefunctions (f> vary continuously, the connection 
factor takes the form 



Tfe = (0(Rfc) I ^(Rfe))exp{ze(i9t(/.(Rfe) | V(Rfc))} + ©(e^) 



(20) 



where the time derivative acts on the wavefunction (f). 

Thus, on multiplying the matrix elements and correction factors together, the following path integral expression 
for the electronic contribution Gna can be obtained as below: 

(*°;Q' I Gna I 4'^;Q> = J I?r(R(T)),7^(R(r))]i?[r(R(T)),0(R(r))](*";Q' I m')) 



(0(Q) I */;Q)exp 



(0(RO \idr~He\ ^{Rr)) 

(0(R,) I i;iRr)) 



dr 



(21) 



The quantities D[il;* = Hjli dfi[tpj{R-i)] are the measures for the coherent state path integral. The time derivative 
is symmetrized by integrating by parts and ignoring boundary terms. Thus, the full path integral between two states 
Q') I 4''') and | Q) | is given by a total action which of the form: 



/ ^1 f'sT (V, o^ , {m)\^^dt -ffe |^(R)) ,^, 
5dR,V^,</>]^y^ {L„(R,R) + ^^(R) I ^(R)^ }dt 

Here the time derivative, 9 t, is the symmetrized derivative defined as: 

(0(R) I zd! I ^(R)) = '-[{dtm) I ^(R)) - im) I dMR))] 



(22) 



(23) 



The purely nuclear contribution to the action integral is contained in the Lagrangian i„ which is defined as follows: 



L„(R,R) = ^MR^-U{R) 



(24) 



Here the quantities R correspond to the set of nuclear coordinates, and M is the mass of the nuclei. U{R) is the bare 
nuclear potential as defined earlier in Eq.(2). The path integral in terms of this action can be written as 

(*°;Q' I Gna I *^;Q) - J i5(^*(R(r),V'(R(r))i?(0*(R(r)),0(R(r))(vl/O | i,{Q')) 

(</.(Q) I VI//) exp [z5t[R(t),V^,0]" 

(25) 

Since, in general, off diagonal matrix elements of the time evolution operator in the coherent state basis are being 
considered, the boundary conditions on the path integral are that the coherent states defined at the end points Q' 
and Q are allowed to vary without constraint. However, in the semiclassical limit, which corresponds to choosing only 
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the subset of continuous paths in coherent state space, the coherent states at the boundaries are constrained to have 
unit overlaps with the initial and final electronic wavefunctions respectively. 

As noted earlier, the necessity of choosing two sets of coherent states at each time slice in the path integral means 
that the paths in general are discontinuous. This is because, in the limit where the discrete time slicing in the 
path integral is made infinitcsimally small, the choice of different states at each such time slicing corresponds to a 
discontinuous jump in the space of coherent states. 

However, the stationary phase limit of the path integral involves continuous paths in state space. Hence, it is 
reasonable to reduce the unconstrained sum over all (discontinuous and continuous) paths to one where the coherent 
states are constrained to be similar to each other. A simple way of doing this is to assume the overlap, (<^(R) | ^p(R)), 
to be nearly unity. This is a reasonable approximation also because contributions to the path integral when this 
overlap is small are negligible. This can be seen from the discrete version of the path integral, Eq.(13) and Eq.(16), 
from which it is clear that as the overlap goes to zero, the exponent in the path integral becomes highly oscillatory, 
and in addition is multiplied by the small overlap. 

The action when the overlap is assumed to be unity can be considered to be a zeroth order approximation to the 
full action, when expanded in terms of the overlap. With this simplification, the path integral action becomes: 

5t[R, V-, = / [in(R, R) + im I ^'dr - H, I ^(R))]dT (26) 
Jo 

This action has a simple interpretation. The time derivative part of the action corresponds to the nonabelian gen- 
eralization of Berry's phase, while the other parts correspond to the dynamical evolution of the system. It is now 
possible to show that the generalized Berry's phase in this action corresponds to nonadiabatic coupling terms between 
Born-Oppcnheimer energy levels. To show this, the Slater determinant states (/> and ip are expanded in a local basis 
of Born Oppenheimer electronic eigenstates as below: 

I V(R)) = $^«^;.[V',R] l*^;R) (27) 

|</.(R))=^^z;^[</.,R] |*^;R) (28) 
ft 

The states | ^^;R) are eigenstates of the electron Hamiltonian 

[ffe(r,R)] I *;.;R) = E^{R) I *^;R) (29) 
In this basis, the phase of the connection factor over the entire path, F can be rewritten as 

r[CR]=^ t dT^{w;{^^;R\idr[w,\^,;K)]-idr[w;{^^;K\]w,\^,;R)) (30) 

This connection factor is the phase obtained by multiplying the connection factors F/j, in Eq.(20) over the entire path. 
A little algebra can be used to demonstrate that the connection factor can be further simplified into the following 
form: 

nCn]= [ dTy^hw;[^,-R]iw^[i;,-R]+h.c)- f dR ■ V A^,(R)«;;[0, R]«;,[^, R] (31) 

Here, the coupling term A^i, is the familiar nonadiabatic coupling vector in the Born Oppenheimer basis: 

A^. = I Vr*,) (32) 
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The second term that contributes to the connection factor has the form of an integral over the gauge potential 
Ap^. In contrast to the case of a simple Berry phase, this "gauge" potential is an infinite dimensional matrix and 
corresponds to a nonabelian geometric phase. Thus the nonadiabatic coupling between Born Oppenheimer surfaces 
has a purely geometric interpretation as arising due to the parallel transport of a local basis of electronic wavefunctions 
along a nuclear path. It is therefore a natural generalization of the Berry phase for parameter dependent quantum 
systemsii. This generalized geometric phase also depends on the off-diagonal elements of the electronic density 
matrix, and consequently couples the slow, nuclear and fast, electronic degrees of freedom. In the limit where an 
infinite number of adiabatic states corresponding to the fast modes can be accessed, this phase represents a dissipative 
coupling between the system of slow nuclear modes and the bath of fast electronic modes^S. 

It should be noted that generalizations of Berry's phase to the case of nondegenerate electronic levels have previously 
been obtained by Mead^^ and by Zygelman i4,2i^ related path integral description of the nonadiabatic dynamics 
has also previously been derived by Moody, Shapere and Wilczeki^. In their work, a path integral was obtained by 
introducing Born-Oppenheimer electronic eigenstates at each nuclear position, in the expression for G„a[R.] given by 
Eq.(12). The electronic matrix elements of GnalR] between initial and final Born Oppenheimer eigenstates | Q') 
and I vPn; Q) are given by 

(*™;Q' I G„a[R(<)] I «'«;Q) = J2 exp[-ie^,„(Q')](^',„; Q' | ^'„,;Ri)exp[-ie^,„,(Ri)] 

mi ,m2 ,. ■ 

(*™i;Ri I *™2;Ri)....(^'„„;RAr | 'J'„;Q) 

(33) 

This expression contains connection factors of the form: 

= (*™,;R, I vI',„^_^,;R,+i) (34) 

It is easy to see that this connection factor can be manipulated as with the coherent state case to give 

= [<5,„„ - <5Rfe • {^,n,;Rk I VR,^,n,+,;Rk)] (35) 

with (5Rfc = Rjt+i — Rfe. This expression can be rewritten in terms of the nonadiabatic coupling to give 

= {eM-^SRk ■ A,„„™,^, (Rfe)]}™„™,+, (36) 

This connection factor is similar to the connection factor, Eq.(20) derived in the coherent state electronic basis. It 
is however harder to handle for the purpose of approximations and computations, due to the discrete nature of the 
summation over Born-Oppenheimer eigenstates. Although the approach presented here is physically equivalent to 
that of Moody, Shapere and Wilczek, the use of electronic coherent states enables semiclassical approximations to the 
nonadiabatic quantum dynamics. This is possible because the restricted sum over only Born Oppenheimer eigenstates 
in Eq.(33) is replaced with an integral over continuous paths in the coherent state space. This provides a practical 
advantage, in that computationally tractable semiclassical approximations to the full path integral become possible. 
It is shown later in this work that stationary phase approximations to the nonadiabatic path integral so obtained, are 
closely related to the Ehrenfest formulation of mixed quantum-classical dynamics. 

The path integral representation derived here is an exact reformulation of the quantum dynamics of the system. 
However, due to the continuum nature of the coherent state basis, it is difficult to numerically evaluate the path 
integral propagator and obtain the full quantum dynamics. 
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III. Analysis 



The path integral derived in the previous seetion ean be used to develop a semiclassical formulation of the nonadia- 
batic dynamics. In this section, stationary phase approximations to the path integral are derived and a semiclassical 
quantization rule is described for the nonadiabatic dynamics. The stationary phase equations so obtained form a 
natural quasiclassical limit for the full quantum dynamics. For the coherent state path integral, this approximation is 
shown to yield the Ehrenfest equations of motion with an additional correction term. This correction term corresponds 
to the geometric contribution to the equations of motion, and is non-zero in general when the nuclear subspace is 
three dimensional or higher. 

For the sake of completeness, the stationary phase approximation to the coherent state path integral is compared to 
a generalized stationary phase type approximation to the nonadiabatic path integral when written out in the basis of 
electronic Born-Oppenheimer eigenstates. It is shown that a generalization of the stationary phase approximation in 
this case yields a quasiclassical theory which is the target of the fewest switches trajectory surface hopping method. A 
natural semiclassical extension of this approximation yields a generalization of the fewest switches method which can 
account for hops in classically forbidden regions. Finally, the derivation of general mixed quantum classical methods 
that combine the Ehrenfest and surface hopping methods is discussed. Thus, this section clarifies and lays down the 
neccesary theoretical foundations for the systematic improvement of mixed quantum classical methods. 



A. Quasiclassical equations of motion 



A stationary phase approximation to the path integral can be made to yield quasiclassical equations of motion for 
the coupled dynamics of the nuclear and electronic degrees of freedom. Such approximations are derived here to obtain 
quasiclassical equations of motion for the system. It is shown that the equations so obtained are a generalization of 
Ehrenfest's equations of motion. 

A density matrix for the electronic subsystem can be defined as 

CT„/3[V', (/>; R] = w*^[(l), IL]wc[iIj, R] (37) 

where the quantities Wa,Wj3 are as defined in Eq.(27) and Eq.(28). To express the nonadiabatic action in terms of 
this density matrix, the explicit form of the action from Eqs.(25) and (30) is written to be 



St[R,(t] = dr i:„(R,R) + ^['w*^'w^ - w*^w^] - R - V A^^u;^[^,t]u;*[^,t] 
Jo ^ ^ 

With this form of the action, the stationary phase approximation corresponds to the condition 



(38) 



JS't[R,(t]=0. (39) 

A detailed derivation of the stationary phase equations is given in Appendix A. at the end of this work. On performing 
the appropriate manipulations, the stationary phase condition for the electronic density matrix, a, can be shown to 
be the following equation of motion: 

i&a/} = (£'a(R) - Eis{R))(7c,i3 + R • ^ [Aa^C7^/3 - (Ta^Ajf}] (40) 
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Here, the energies S„(R), -B/3(R) are the eigenvalues of the Born Oppenheimer eigenstates as defined in Eq.(29). The 
corresponding equation of motion for the nuclear coordinates, R with mass M and the bare nuclear potential U (R) 
can be written to be 

MR = -Vr?7(R) - J2 Vr^^(R)ct^,. + Yl A^^'^M. - '^''-^ ^ ^ ^-'*) (^1) 

The stationary phase equation for the nuclear coordinates, R, contains a term that is proportional to the time 
derivative of the off diagonal electron density matrix cr^^. In addition to the other terms which are proportional 
to the off-diagonal density matrix, this is an explicit contribution to the force acting on the nuclei due to the time 
variation of the electronic coherence. Hence, this accounts for changes in the effective nuclear force due to electronic 
dephasing. 

Eq.(41) can be further simplified and written explicitly in terms of the elements of the electronic density matrix, a 
and the nonadiabatic couplings A as, 

MR = -VnV - Y "^nE^a^f, + i ^ E^^A^^a^^ - ^ cr^^R x (Vr x A^^) + 

If a nonabelian "vector potential" is defined as an infinite dimensional matrix with elements 

[A(R)]^, = A^, (43) 

then the corresponding gauge invariant, nonabelian "magnetic field" is given by 

B^^ = Vr x A^^ - i(A X A)^^ (44) 

Prom these definitions and Eq.(42), the stationary phase equations for the nuclear coordinates can be written in the 
form 

-E^!^ (R)A^^ cTy^ - R X ^ Byjucr^^ (45) 

Eq.(40) corresponds to the familiar Liouville equation for the electronic density matrix, while Eq.(45) is diff'erent 
because it contains corrections due to the nonadiabatic coupling. The first three force terms correspond to the 
eff'ective force in the Ehrenfest formulation of mixed quantum classical dynamics, while the last term is a geometric 
correction to the Ehrenfest equations. This term therefore constitutes a new addition to the equations of motion and 
will need to be retained for consistency with the full quantum dynamics. When this term is ignored, the equations of 
motion so obtained correspond to the Ehrenfest equations of motion. 

The equation, Eq.(45), for the nuclear coordinates satisfies the invariance under general unitary transformations of 
the electronic basis. To sec this, note that the nonadiabatic coupling matrix A(R) is defined through the nonadiabatic 
coupling matrix elements as 

[A(R)]^, = A^,(R) = I Vr*.) (46) 

To write this in a more transparent form, infinite dimensional column vectors * can be constructed with elements 
given by the electronic basis wavefunctions The matrix can then be rewritten as a product of column vectors 

A(R) = *t(_^VR)* (47) 
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It is straightforward to show that under a unitary transformation, n(R) of the electronic wavefunctions, 



* =n(R)* 



(48) 



the nonadiabatic coupUng matrix transforms as 



A(R) ^ nt An - mtvRn 



(49) 



The "magnetic field", B with matrix elements B^^ correspondingly transforms as 



(50) 



The density matrix elements a^^^ also transform similarly under the unitary transformation, and hence the " magnetic 
field" contributions to the equations of motion are invariant under such transformations. This invariance is a general 
version of the the usual invariance under gauge transformations, of vector potentials and magnetic fields in classical 
electromagnetism. It should be noted that the unitary transformations under which the electronic wavefunction basis 
transforms are infinite dimensional matrices, and hence the group, U{oo) that these unitary transformations belong 
to is the group of symmetries for the nonadiabatic corrections. On the other hand, the terms in the equations of 
motion which are gradients of the electronic energy eigenstates break this symmetry, since the energy eigenvalues 
are, in general, not completely degenerate. When degeneracies exist for a subset of the electronic energy levels, the 
contribution to the nuclear force from the degenerate states is symmetric under unitary transformations corresponding 
to the degenerate subspace. The group corresponding to this symmetry is U (n), where n is the degree of the degeneracy. 
When this symmetry group exists, the time evolution of the degenerate electronic subspace picks up a Berry phase^^ 
which contributes to the nuclear force in Eq.(45) in a fashion identical to that of the nonadiabatic coupling between 
nondegenerate states. Thus, the equation of motion for the nuclei, Eq.(45), is valid irrespective of whether the states 
in the electronic subspace are degenerate or not. For example, in the case of a conical intersection between two 
potential energy surfaces, the wavefunction becomes multivalued. To resolve this multivaluedness, a geometric phase 
due to a vector potential can be added to the wavefunctions. The geometric forces in Eq.(45) are the classical analog 
of this procedure. 

It should be noted that the geometrical, "magnetic field" contributions in Eq.(45) are zero for one dimensional 
nuclear potential energy surfaces. This can be intuitively understood through the discussion in Sec. II. A. This magnetic 
contribution is due to the "rotation" of the local electronic Hilbert space along the nuclear path. For one dimensional 
nuclear potential energy surfaces, this "rotation" does not have any torsional components, and is simply a pure 
parallel transport of the basis. Thus, for single nuclear dimensions, these forces are zero and do not contribute to 
the dynamics. However, for potential energy surfaces parametrized by higher dimensional nuclear coordinates, these 
forces are non-zero and of importance to the nonadiabatic dynamics. These couplings therefore cause an effective 
magnetic force contribution to be added to the usual Ehrenfest force, in the mixed quantum classical dynamics. The 
magnetic forces are the quasiclassical analog of the additional Berry-Mead phase which the electronic wavefunctions 
are multiplied by, to maintain their single valuedness. 



A complementary description of the full nonadiabatic quantum propagation of the nuclei can be obtained by writing 
the nonadiabatic propagation operator G'„a[R] from Eq.(12) in a basis of electronic Born Oppenheimer eigenstates. 



B. The Discrete Path Integral 
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As mentioned previously such a path integral expression for the propagator has been written down by Moody, Shapere 
and Wilczek. Explicitly, if the wavefunction for the combined system is given by 

(Q|*,t) = 5^^„(Q,i) |*„;Q) (51) 
then the time evolution of the nuclear states can be written as 

UiQ', t) = Y.j KmniQ!, t I Q, 0)e„(Q, 0)dQ (52) 

n 

The propagator Kmn is the matrix element of the time evolution operator when evaluated between an initial electronic 
Born-Oppenheimer state, defined for an initial position Q and a final electronic Born-Oppenheimer state, 4"^, 
defined for a final nuclear position Q'. This matrix element can be evaluated using the Trotter product formula to 
slice the time t into infinitesimal time intervals, and introducing a basis of electronic Born Oppenheimer states at 
each time slice. The resulting propagator is the same as the propagator defined in Eq.(33). 
Explicitly, this propagator is the matrix element 

K„^n(q,\ i I Q, 0) = Q' I e-'"' \ Q) (53) 
If the Born-Oppenheimer limit is assumed to be valid, then the propagator is diagonal in the indices m, n. 

Ktl{Q',t I Q,0) = 6mnKt"m{Q',t | Q, 0) (54) 
Here, if^^ is the nuclear propagator in a potential given by Kn(Q) = U{Q) + Em{Q) 

KmmiQ', i I Q, 0) « (Q' I e-*(^"+^-(^))* I Q) (55) 

The nonadiabatic nuclear propagation can be thought of as a sequence of propagation along Born-Oppenheimer paths 
with the propagator defined above in Eq.(54), and "hops" occurring in between to switch between Born-Oppenheimer 
paths. This is illustrated in Fig. 2. Detailed considerations of the structure of the paths that contribute to the 
transition amplitude Kmn thus point to the origins of surface hopping type approximations. 

To derive explicitly quasiclassical techniques based on the discrete path integral, consider the effective propagation 
operator defined in Eq.(12), which corresponds to matrix elements of the full propagation operator between initial 
and final nuclear position eigenstates. To recapitulate, the matrix elements are given as an integral of a path ordered 
exponential Gna [^{t)] ■ For a given nuclear path R(f) the matrix elements of Gna between an initial Born-Oppenheimer 
electronic state and a final state are given by 

Q' I G„„[R(i)] I Q; *n) = (Q'; *m I e-^«(Q')^e-^«(^^)^..e-^«(Q)^ I Q; (56) 

Prom this relationship, the quantum nonadiabatic propagator between an initial Born oppenheimer electronic wave- 
function I ^'„; Q) defined at an initial nuclear position Q and a final Born Oppenheimer electronic wavefunction, 
I ^m', Q') defined for a position Q' can be written as 

Kmn{Q',t I Q,0) = J Dn{T){Q';^^m \ G„a[R(t)] | Q; (57) 

If a nonadiabatic transition takes place during the interval {t', t' + e) for sufficiently small values of e, the propagator 
has the composition property 

KmniQ!, i I Q, 0) = / rfR'rfR" KmpiQ', t \ R', t'){^p; R' | R")^s:g„(R", i' + e | Q, 0) (58) 

pq 
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The overlap between electronic wavefunctions at R' and R" is the connection factor, rpq(R',R"), discussed in the 
previous section, and is responsible for the nonadiabatic coupling between levels. As before, it can be rewritten as: 

{■^p; R' I 4-,; R") = rp,(R', t \ R", t') « (<5p, - zA^, • (R' - R")) (59) 

with the nonadiabatic coupling A,„„ given by Eq.(32) as before. This composition law for the exact nonadiabatic 
propagator can be used to derive a Born series expansion for the propagator in terms of the Born Oppenheimer 
propagator. The physics behind this is illustrated in Fig. 2. 

K,nn{Q',t\ Q,0) = K'mUQ',* \Q,0) + J dR'K^UQ'^t I R',i')r™„(R',R")<°(R",i' + e I Q,0) 

+terms with two nonadiabatic transitions + .... 

(60) 

Here, the propagator is written as a sum over paths where there are successively higher number of nonadiabatic 
transitions. To obtain a classical nuclear propagation scheme, the nuclei can be propagated along the stationary 
phase solutions for those pieces of the nuclear path where the propagator is of the Born Oppenheimer type. In 
other words, the Born- Oppenheimer limit of the nonadiabatic propagator can be approximated with a semiclassical 
expansion around the stationary phase path as discussed below. 
It is easy to see that along a Born Oppenheimer path, 

.0^ = i?,„^m (61) 
ot 

Furthermore, if a single nonadiabatic transition occurs from n ^ m on a sufficiently short timcscale e, the amplitude 
S,m changes by 

<5ern(R',t) - -z^<5R' • A„„(R')e„(R',0 (62) 

n 

Putting the two equations together, the total time derivative of the wavefunction ^ is recovered to be: 

itn(R, t) = i/„^„ + R- • I] A„p(R)ep(R, t) (63) 
p 

It is important to note that this equation is valid under the assumption that nonadiabatic transitions occur on a very 
short timescale in comparison to the timescales that govern the rest of the dynamics. In general, there is a nonzero 
probability for nonadiabatic transitions occuring over longer timescales and accompanied by finite changes, JR' in the 
nuclear coordinates. This equation is derived assuming that the contribution from such transitions is small, and that 
dominant features of the nonadiabatic dynamics can be captured by considering transitions that occur on infinitesimal 
timescales. 

To derive a semiclassical limit for the nuclear evolution, the nuclear wavefunctions ^„(R, f) are approximated as 
Gaussian wavefunctions centered around the instantaneous classical nuclear path on the Born Oppenheimer surface, 
according to the prescription of Heller—. 

C„(R, t) « c„(t)i^„[R„(t), P„(t)] (64) 
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The functions Fn(t) are instantaneous Gaussian wavepackets centered around the classical Born-Oppenheimer trajec- 
tories given by 

Furthermore, this Gaussian wavepacket approximately satisfies (when the Born Oppenheimer potential energy, U (R)-l- 
En(R), is expanded to quadratic order in R) the time dependent Schrodinger equation: 

zF„ « i/„F„ (66) 

Substituting the ansatz, Eq.(64), into Eq.(63), one obtains for the time evolution of the coefficients, c„ 

p 

Making use of the approximation in Eq.(66), the time evolution of c„ becomes, in the semiclassical limit: 

ic„Fn = X! ^ ' ^npFpCp (68) 

p 

Multiplying by F* and integrating over the coordinates R gives 

icn - ^ (^^„ I R • A„p I Fp)cp (69) 
p 

This can be rewritten in terms of a density matrix cr„m = c„c^ as 

i^nm — ^ ^ [(-^n | R ' ^np \ Fp)<7p^ f^np{Fp \ R ■ -^pra \ -^m)] ('^0) 
P 

Furthermore, if the Gaussians Fn are assumed to be delta functions in the positions R, the above equations simplifies 
further into: 

^^nm ~ ^ ^ [R ' ^np^pni f^npR ' -^pm] ('^-^) 
P 

where the variable R and its velocity, R are classical quantities corresponding to the center of the Gaussian function 
Fn- Since the functions F„ are Gaussians centered on the classical nuclear Born-Oppenheimer trajectories, the classical 
variables R satisfy Newton's laws: 

MR = -Vr[C/(R) + £;„(R)] (72) 

when the trajectory being considered is that belonging to the n*'' Born-Oppenheimer surface. Further, the rate of 
change of the diagonal elements, Gnn is given by 

i^nn R ' ^ ^ \-^np^pn ^np^pri\ ('^*^) 
P 

It is clear from this equation and the fact that the nuclear propagation always occurs on a single Born-Oppenheimer 
surface, that these equations describe the surface hopping method-'^ , for which the fewest switches algorithm is a 
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Monte Carlo solution. Thus, the fewest switches surface hopping (FSSH) approach^ corresponds to a semiclassical 
approximation to the exact quantum nonadiabatic dynamics, similar to the Ehrenfest method described in Sec III. A. 

It is also interesting to note that Eq.(70) can be used instead of Eq.(73) to generate the surface hopping criterion. 
Eq.(70) provides for a nonzero hopping probability even in regions where transitions are classically forbidden, unlike 
the usual FSSH method. A similar extension has also been proposed previously by Neria and Nitzan^^, on physical 
grounds through the construction of a semiclassical golden rule. In addition to this, it is clear from this discussion that 
semiclassical propagation of the nuclear wavefunction, either using the Gaussian form, Eq.(65) and Eq.(66) proposed 
by Heller—, or by the use of a different basis, lead to generalized schemes for nonadiabatic dynamics that go beyond 
the mixed quantum classical limit. The use of frozen gaussians within this framework correspond to a scheme for 
nonadiabatic dynamics that is similar to the multiple spawning method proposed by Martinez and co-workers^iiS^. 

The approximation of hops occurring over infinitesimal timescales that is required to derive both Eq.(70) and 
Eq.(73) causes phase changes in the wavefunction due to the hops to be neglected. Furthermore, hops need not 
occur on infinitesimal timescales, and the quantum nuclear paths which contribute to the dynamics are in general 
not differentiable, rendering the nuclear velocity to be a poorly defined quantity. In circumstances where such paths 
contribute significantly to the dynamics surface hopping type methods will be invalid. Additionally, in the case where 
two (or more) electronic potential energy surfaces become degenerate, the nonadiabatic coupling also picks up nonzero 
diagonal elements which could act on the dynamics along significant sections of the nuclear path. In such situations, 
an approximate means of treating the dynamics would be to add the diagonal terms in the nonadiabatic coupling to 
the nuclear Hamiltonian and perform the nuclear propagation on each Born Oppenheimer surface with an additional 
magnetic force due to the diagonal terms in the nonadiabatic coupling. The off diagonal nonadiabatic couplings can 
then still be treated as acting instantaneously to cause hops. It is likely that such modifications could be of value in 
improving the accuracy of the FSSH method in such situations. 

Finally, it is useful to discuss the relative utility of the discrete basis path integral presented here versus the coherent 
state approach that leads to the Ehrenfest equations in the previous section. The natural stationary phase limit of the 
coherent state path integral is the mean field Ehrenfest theory, while the discrete path integral has a corresponding 
stochastic limit. These two limits can be thought of as corresponding to two different types of statistics in the 
hopping events. The coherent state path integral is a better starting point for approximations when there is significant 
nonadiabatic coupling between multiple Born-Oppenheimer surfaces, i.e when hops can occur with high frequency. In 
such a situation, a mean field treatment is a good first approximation to the dynamics. However, the discrete basis 
approach is better suited to studying the dynamics in the presence of nonadiabatic dynamics between a few levels, and 
when Born Oppenheimer surfaces begin to diverge. In such a situation, the hopping events are relatively infrequent, 
and nonadiabatic behaviour is sensitive to the detailed structure of neighbouring Born-Oppenheimer surfaces. The 
discrete basis approach, for which the surface hopping algorithm is an approximation is a good starting point for 
approximation in this situation. 

C. Semiclassical quantization conditions 

The nonadiabatic action Eq.(38) can be used to describe a semiclassical quantization condition when the system is 
in a bound state. This condition is derived in a fashion similar to that used to study the semiclassical quantization 
of coherent state path integral a^^'^^ , and in quantum field theories^^. 
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If an assumption is made that the system executes periodic orbits, then a quantization condition can be written. 
To derive this quantization condition, the trace of the propagator is written as 

K„,{T)=Tr[e^p{-iHT)] (74) 

The trace is related to the density of states as below. 

If the trace corresponding to p{E) is expanded in terms of the eigenstates of the entire (the combined nuclear and 
electronic) system, p{E) can be written in terms of the system energy eigenvalues En as 

p(E) = Y (76) 

> E-En + iti ^ ' 

It is clear from this that the poles of the function p{E) correspond to the energy eigenvalues of the system. Hence, 
the poles of the Fourier transform of the propagator trace are related to the eigenstates of the combined system in 
this fashion. 

The Fourier transform of the trace of the propagator can be written in terms of the nonadiabatic coherent state 
path integral as follows: 

p{E) =ij^ dT J dUd-^ J D[Il{t)]D[ip{t)] exp {iSrlB,, ip, cj)] + ET) (77) 

Here, the quantities ■(/;,(/) are coherent state basis elements and the action St is the same as defined in Eq.(26). To 
obtain an approximate quantization rule, only diagonal coherent state matrix elements in Eq.(26) are considered. 
Eq.(75) which relates the Fourier transform p{E) to the energy spectrum of the system also implies that p{E) is 
singular when the energy E is equal to any of the system eigenenergies. This implies that the path integral of the 
propagator in Eq.(77) is a maximum for such values of the energy E. The propagator also acquires its maximal values 
when its phase is stationary. Consequently, the density of states, p{E), is dominated by contributions to the path 
integral from stationary phase paths. 

Furthermore, if it is assumed that periodic orbits exist for the combined dynamics of the system, then the density 
of states in Eq.(75) can be approximately written as 

p[E] = ^ Cfe exp [iS[n{E)] + iETk{E)] (78) 

p.o 

where the time T{E) is the time period for the classical traversal of the fcth periodic orbit. Since, multiple traversals 
of a given periodic orbit also constitutes a periodic orbit, this sum can be simplified into a sum over periodic orbits 
p, which cannot be reduced further into smaller sets of periodic orbits: 

r,lF^-S- ,r ^M^S[T,{E)]+^ET,{E)] 

This expression has poles given by 



27r(n +\) = ET{E) + 5[T, E\ (80) 



This quantization condition can be explicitly written out to be 

i-T, 



ETr.. 



{E)+ f " [i„(R, R) + (V- I idt - He I V)l = 27r(n + J) (81) 
Jo ^ J i 
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It should be noted that this quantization condition is vahd under the assumption that periodic orbits corresponding to 
the effective classical system defined by the action S exist. For nuclear motion coupled to a large number of electronic 
levels, finding the periodic orbits of the system is a difficult task. However, for small values of the nonadiabatic 
coupling, this quantization rule provides an approximate means of estimating the total energy spectrum of the system, 
given an initial set of Born-Oppenheimer surfaces. Furthermore, if the nonadiabatic transitions couple only a finite 
set of electronic Born-Oppenheimer surfaces, this rule may have validity because periodic orbits for the combined 
nuclear and electronic density matrix variables could exist. 

This rule can be written in terms of the total connection factors r[CR] for the closed path Cr as 



A special case of the Born Oppenheimer problem, for which the semiclassical quantization rule is useful, is that of 
a single vibrational coordinate, p, q, coupled to a set of electronic levels. In this situation, the quantization rule reads 



In the absence of the nonadiabatic couplings, A„,„, the quantization rule gives the exact vibrational energy levels 
(assuming that the vibrations are purely harmonic). The contribution due to nonzero nonadiabatic coupling can be 
interpreted as a shifting of the energy levels of the oscillator to account for energy exchange with the electronic bath. 
If an orbit is executed, nonadiabatic transitions between vibrational states can occur, which leads to energy transfer 
between the oscillator and the electronic system. 



This work has presented an exact first principles path integral formulation of quantum nonadiabatic dynamics. The 
path integral approach developed here has the advantage that it enables computationally tractable stationary phase 
approximations to the full path integral. Previous path integral based approaches to this problem have stationary 
phase approximations that lead to nonlocal theories of mixed quantum classical dynamics^iS., or require considering a 
constrained set of paths for the purposes of approximation^^'^^. The method presented here therefore represents an 
advance over previous approaches in this regard. 

A second consequence of this work is that corrections to the quasiclassical nuclear dynamics have been obtained. 
These corrections reflect underlying geometric features of the electronic Hilbert space. The corrections so obtained 
have the advantage that they can be calculated through straightforward extensions of current mixed quantum-classical 
simulation methodologies. Applications of this approach to study the dynamics of real systems in the presence of 
conical intersections and other degeneracies in the electronic energy levels are currently being considered. Finally, 
semiclassical quantization rules have been proposed for nonadiabatic dynamics in the quasiclassical limit. These rules 
are potentially of utility in determining approximate energy spectra when the Born-Oppenheimer approximation is 
not valid. 

A third consequence of this work has been to provide a link between exact quantum theories of nonadiabatic dy- 
namics and semiempirical surface hopping methods. This work establishes a rigorous theoretical basis for the surface 
hopping method. It has been demonstrated that the fewest switches surface hopping approach can be regarded as a 




(82) 




k,l 



(83) 



IV. Conclusions 
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generalized stationary phase approximation to the exact quantum dynamics, in complement to the Ehrenfest formu- 
lation of mixed quantum classical dynamics. Natural extensions of the surface hopping method to include classically 
forbidden nonadiabatic transitions, as well as the inclusion of nuclear quantum effects through the construction of 
rigorous multiple spawning type algorithms have been demonstrated. This work also allows for the construction of 
semiclassical golden rules to describe nonadiabatic transitions. 

The present work also provides a clear perspective on issues concerning equilibrium between the quasiclassical 
nuclear and quantum, electronic subsystems. It is clear from the theory presented here that traditional Ehrenfest 
methods do not satisfy equilibrium conditions if only single mean field trajectories are considered. Each such trajectory 
corresponds to a stationary phase solution for the nonadiabatic propagator at a given set of initial conditions. How- 
ever, thermal properties of the system enter through sampling of initial conditions. Hence, if appropriately sampled 
ensembles of mean field trajectories are chosen, a mixed quantum-classical equilibrium can be established. Correct 
sampling of mean field trajectories will require an initial sampling of the reduced density matrix for the nuclei and 
should also correctly account for initial phase correlations according to Eq.(26). It is possible that observed violations 
of mixed quantum classical equilibrium by the Ehrenfest method are due to an incorrect sampling of cither of these 
contributions. Since the surface hopping approximation is obtained by looking at the propagation of the system in 
each of the states occupied by it, an ensemble averaging of surface hopping trajectories will sample the initial density 
matrix of the system, and hence this approach and any systematic generalizations based on it can be expected to 
reproduce the correct mixed quantum classical equilibrium. This has been demonstrated through a separate set of 
arguments^!. 

This work provides a different conceptual perspective on the nature of nonadiabatic dynamics. Nonadiabatic 
couplings can be regarded as geometrical objects which are generalizations of Berry's geometric phase. Thus, they 
can be regarded on the same footing as the geometrical phases that arise due to conical intersection o^^'^^'^° . This 
similarity between nonadiabatic couplings and geometric phases that arise in degenerate systems has also been noted 
by previous author a^^'^^ . In principle, the methods described in this work do not differentiate between systems which 
have degeneracies in energy and those which are nondegenerate. Therefore, the nonadiabatic path integral formulation 
presented here is of considerable generality and may be of practical utility in computational studies of nonadiabatic 
dynamics. 

A potentially useful approach to study exact quantum nonadiabatic dynamics is to consider quantum Monte Carlo 
schemes to evaluate the coherent state path integral. Direct applications of quantum Monte Carlo schemes are likely 
to face difficulties due to the sign problem. However, simplifications could be made by using an effective single 
electron ansatz to replace the matrix elements in the action Eq.(26). Such an ansatz could be rigorously constructed 
by using the Runge-Gross theorem'^^ which relates the exact one electron time dependent density to the expectation 
value of the electronic Hamiltonian in Eq.(26). Extensions of the time dependent density functional approach to one 
electron density matrices can then be used to accurately evaluate the path integral due to the action in Eq.(26). This 
has the added advantage of naturally incorporating time dependent density functional theory with the nonadiabatic 
propagation of a quantum system. This approach will be explored in later work. 

Although the exact coherent state path integral presented here is a challenge for numerical computation, it is 
possible to develop numerically more tractable semiclassical propagators that describe the nonadiabatic dynamics as 
an approximation to the coherent state path integral. One such approach would involve a stationary phase evaluation 
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of the electronic dynamics (the Ehrenfest equations for the electron density matrix) , while propagating the nuclear 
motion using a semiclassical initial value representation. More promising applications of this work are in developing 
better approximate methods to describe quantum nonadiabatic dynamics. One potential application is to develop 
methods which systematically combine Ehrenfest and surface hopping type algorithms. To enable this the propagation 
of the nuclear dynamics can be divided into sections wherein different basis sets are used depending on physical 
considerations regarding the nature of nonadiabatic interactions along the nuclear paths. Another natural extension 
is in refining the surface hopping method to include classically forbidden nonadiabatic transitions by approximating 
the nuclear wavefunction with a gaussian wavepacket. Alternatively, more general semiclassical propagation schemes 
can be used to propagate nuclear dynamics in the discrete representation of the path integral. This would facilitate 
the inclusion of nuclear quantum effects as well as more accurate treatment of dynamics in the presence of strong 
nonadiabatic coupling between Born-Oppenheimer surfaces. 

The use of coherent states, and the derivation of the path integral presented here is analogous to path integral 
treatments of quantum spin systems. This analogy with quantum spin systems also forms the basis for classical 
mapping techniques^'^i^'^i^^ that have been developed to study mixed quantum classical dynamics. Finally, the 
path integral methods presented here are apparently amenable to perturbative approximations around its stationary 
phase value. However difficulties remain due to the nature of the nonadiabatic corrections to Born Oppenheimer 
dynamics. The nonadiabatic coupling between Born Oppenheimer levels consists of matrix elements of the nuclear 
momentum operator. The nuclear momentum operator is unbounded in the Hilbert space, and further, the nuclear 
paths that contribute to the path integral for the system are in general not differentiable. Hence, a perturbation 
expansion in the nonadiabatic coupling element is not well defined in the path integral framework^^. In the case of 
Ehrenfest dynamics, this issue also manifests through the assumptions made in constructing the coherent state path 
integral. The assumption made in coherent state path integrals is usually that the paths chosen in the coherent state 
space are continuous. This is not valid in general for the problem of nonadiabatic dynamics. Similarly, extensions of 
surface hopping type methods will need to account for non-differentiable nuclear paths where the nuclear velocity is 
not defined. A means of overcoming this problem for the coherent state path integral is to regularize the path integral 
by considering higher order time derivatives in the action, in analogy to the procedure followed for spin system a^^'^^ . 
A perturbative treatment of nonadiabatic dynamics that takes these issues into account is being worked on and will 
be presented at a later date. 
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Appendix A. Coherent state path integral: Stationary phcise equations. 

The derivation of the stationary phase equations in Sec.III A. is described here. The effective action in the coherent 
state representation is given by 

St^ f LrdT (Al) 
Jo 

with a Lagrangian Lt given by 

^ MR2 i 

Lt = Y^ - [/(R) - ^ E„, {Il)w*^W^ + - [w*^Wm - W*^Wrn] " R ' ^rnnW*rnWn (A2) 

7=1 m mn 

For the stationary phase path, the action is stationary when varied with respect to the variables R and Wm, Wm which 
determine the path. In other words the stationary phase condition: 

5St = (A3) 

implies the equations: 

^^-^ = (A4) 

dt 9R aR ^ ' 

^=0;^ = (A5) 

dWm SW^ 

The first of the three equations above is the usual Lagrangian equation of motion for Newtonian mechanics which 
governs the classical dynamics of the variable R(t), while the next two equations render the action stationary with 
respect to variations of the electronic density matrix. The equations for the variables Wn,Wn ^■re obtained by using 
the appropriate unsymmetrized form of their kinetic energy term (the fourth term in the right hand side of Eq.(A2)). 
On applying these equations of motion to the Lagrangian, one obtains 

^{MR - ^ A„„CT„„} = - Vr[/ - ^ o-„j„ VR-Er„(R) - ^ Vr(R • A„„)o-„„ (A6) 

m,n m m,n 

-iwl = R • 5^ <A„„ + En{K)wl (A7) 

m 

= R • ^ AnmWm + -E„(R)w„ (A8) 

m 

The second and third equations, Eq.(A7) and Eq.(A8), can be combined to give the Liouville Von-Neumann equations 
for the electronic density matrix, a^rn = WnWm- 

i^nm = {En — Em)(^nm + R- • ^ [^nkC^km — (^nk^km] (A9) 

k 

Further, the identity 

Vr(C • D) = (C • V)D + (D ■ V)C + C X (V X D) + D X (V X C) (AlO) 
can be used to simplify the last term in Eq.(A6) to 

Vr(R • Anm) = (R • V)A„„ - R X (V X Anm) (All) 
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Using this, the equation of motion for the nuclear coordinates becomes 

MR = ^ [A„„(T„TO + Amn&nm] " Cr^m VR^m(R) - VrJ7(R) 

m,n 

The time derivative of A can be rewritten as Amn = (R- • V)Am„ and hence cancels with the corresponding term on 
the right hand side. Thus, the equation of motion simplifies to: 

MR=-J2 <^mm Vr^™(R) - VrC/ (R) + [AmnCrnm - CT„™R X (V X A„„)] (A13) 

m m,n 

Furthermore, the time derivative of the electronic density matrix, &nm, can be replaced by using the Liouville-Von 
Neumann equations of motion for the electronic density matrix, to obtain 

MR = -VrC/ (R) - ^ £7„„Vr£;„(R) + i ^ E^nAmn<ynm - 

m m,n 

i A„„{R • [a, A]nm} - Y ^rim^ X (Vr X A„„) (A14) 

mn mn 

The first three terms are the standard Ehrenfest force acting on the nuclear variables, while the last two terms are 
new. Further, using the vector identity 

B X (C X D) = (B ■ D)C - (B ■ C)D (A15) 
The penultimate term on the right hand side can be simplified and rewritten as 

iA„„{R • [a, A]nm} = -icrnmR X (A X A)mn (A16) 
putting this together with Eq.(A14) and defining an effective magnetic field by 

B„„ = Vr X Amn - i{A X A)„„ (A17) 
the equation of motion for the nuclear coordinate is obtained as 

MR = - VrC/ (R) - Y '^mm'^REmiR) + i ^ i^mnA„„CT„„ - R X ^ (T„„B„„ (Al8) 

This is the equation of motion discussed in the text of the paper, and is the Ehrenfest equation of motion with an 
additional contribution from the "magnetic field", B. 
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FIG.l: A schematic describing parallel transport of the electronic Hilbert space along a given nuclear path. The orthogonal 
axes represent the local electronic Hilbert space at a given nuclear position. As the nuclear path is traversed, the Hilbert space 
is transported along with it leading to the factor,r, which connects the Hilbert spaces at different nuclear positions Ri,R2 
while R, Q are initial and final positions for the given nuclear path. 

FIG. 2: A pictorial representation of the nonadiabatic paths that contribute to the quantum dynamics. Two different paths 
corresponding to two nonadiabatic transitions are shown. The propagation along the horizontal lines is governed by the Born 

Oppenheimcr propagator, while vertical columns represent nonadiabatic transitions. The widths, ti, on the x-axis of the 
columns are the times over which nonadiabatic transitions take place. In the surface hopping approximation, these widths 
become infinitesimally small. 
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